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Abstract 

We derive an equation to compute directly the expected occupation 
time of the centered Ornstein-Uhlenbeck process. This allows us to iden- 
tify the parameters of the Ornstein-Uhlenbeck process for available oc- 
cupation times via a standard least squares minimization. To test the 
method, we generate occupation times via Monte-Carlo simulations and 
recover the parameters with the above mentioned procedure. 

1 Introduction 

Nonwoven materials or fleece are webs of long flexible fibers that are used for 
composite materials, e.g. filters, as well as in the hygiene and textile industries. 
They are produced in melt-spinning operations: hundreds of individual endless 
fibers are obtained by the continuous extrusion of a molten polymer through 
narrow nozzles that are densely and equidistantly placed in a row at a spinning 
beam. The viscous or viscoelastic fibers are stretched and spun until they 
solidify due to cooling air streams. Before the elastic fibers lay down on a 
moving conveyor belt to form a web, they become entangled and form loops 
due to the highly turbulent air flows. The homogeneity and load capacity of 
the fiber web are the most important textile properties for quality assessment of 
industrial nonwoven fabrics. The optimization and control of the fleece quality 
require modeling and simulation of fiber dynamics and lay-down. Available 
data to judge the quality, at least on the industrial scale, are usually the mass 
per unit area of the fleece. 

A stochastic model for the fiber deposition in the nonwoven production was 
proposed and analyzed in Ref. |BGKMW08l IGKMW07] . Its core is a stochastic 
Ornstein-Uhlenbeck process for the random motion of the fiber. The aim of 
this paper is to determine the parameters of the Ornstein-Uhlenbeck process 
from available mass per unit area data, i.e. the occupation time in mathematical 
terms. For the sake of simplicity, we focus on a one-dimensional version of the 
Ornstein-Uhlenbeck process. 

The paper is organized as follows: In Section [5] we introduce the Ornstein- 
Uhlenbeck process as a prototypic model for the fiber deposition. Section [3] 
is devoted to the derivation of the expectation value for the occupation times. 
An algorithm to estimate the parameters in the Ornstein-Uhlenbeck process 
from available occupation times is presented along with numerical experiments 
in Section H) Finally, we draw some conclusions and give an outlook to open 
questions. 

2 Model 

As a prototypic model for the fiber deposition, we consider the general one- 
dimensional Ornstein-Uhlenbeck process 

dU t = \{n - U t )dt + adW t , U = U(0) el. (1) 
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The equilibrium [i € R and the stiffness A > govern the deterministic part 
whereas the diffusion parameter a > and a standard Wiener process (or 
Brownian motion) Wt contribute the stochastic part. For the sake of simplicity, 
we consider mainly the centered process X t satisfying 

dX t = -XX t dt + adW t , X = . (2) 

The real- valued random variable X t models the deposition point of an indi- 
vidual fiber on the fleece. If we follow the random variable over a time interval 
[0,T] for T > 0, we obtain the path of an individual fiber. To introduce the 
mathematical analog of the mass per unit area we need the following definition. 

Definition 2.1 (Occupation time). Let T > and consider an interval [a,b] C 
R, where a = — oo or b = oo are allowed. The occupation time Mx,[ a ,b] * s 
defined as 




— x) dx dt . 



Here, 1[ ,&] denotes the indicator function of the interval [a,b] and So(X t — x) 
is the Donsker's delta function introduced in Definition \3.8l below. 

Remark 2.2. The occupation time is a random variable itself. It models the 
time, the random process spends inside the spatial interval [a, b] during the time 
interval [0, T] . In terms of our physical model for the nonwoven production, the 
occupation time can be interpreted as the mass of fiber material deposited inside 
the interval [a,b], i.e. the mass per unit area of the final fleece. This quantity 
is easily accessible even on the scale of industrial production and hence it will 
serve as the input to our parameter estimation problem. 

In the next chapter, we will present tools from white noise analysis to derive 
the expectation of the occupation time for the centered Ornstein-Uhlenbeck 
process X t . Although it is possible to derive the results by classical stochastic 
analysis methods, we use a white noise approach to generalize the concepts also 
to higher dimensions, where one can give a rigorous meaning to multidimen- 
sional Donskers Delta functions as a white noise distribution, in later research. 
Moreover in future work an extension to more complicated processes (e.g. with 
fractional noise term) is planed. Thereafter we show, how to estimate the pa- 
rameters A, a of the process from available data for the occupation times. 

3 Theory 

We start by considering the Gel'fand triple 5(R) C i 2 (R) C S"(R), where S(R) 
denotes the Schwartz space of rapidly decreasing smooth functions, L 2 (M.) the 
Hilbert space of real- valued square integrable (equivalence classes of) functions 
on R w.r.t. Lebesgue measure and 5"(R) the topological dual of S'(R), i.e. the 
space of tempered distributions. This particular choice is the usual one in 



3 



white noise analysis [HKPS93 . By (/, w) we denote the duality pairing between 
ui £ S'(R) and / £ 5(R), an extension of the standard inner product on L 2 (R) 
in the sense of a Gel'fand triple. 

Next, we want to introduce a probability measure on the space 5'(R). 
Therefore, we consider the cr-algebra £?(5'(R)) generated by the cylinder sets 
{(/, ■) : / £ 5(R)}. The white noise measure \i on (S'(R),B) is given via Min- 
los' theorem [BK951 IH18OI IHKPS93j by its characteristic function C 

/ exp(»(/, W ))d / iH=expf-i|/| 2 ') =C(f) 

JS'(R) \ 1 J 

for f £ 5(R). 

Remark 3.1. The Hilbert space of complex-valued square-integrable functions 
w.r.t. this measure /j, is denoted by L 2 (/i) = L 2 (5'(R), B, /i). For f,g £ 5(R) we 
have the isometry 

/ (9,w) dfi(uj) = / f(s)g(s)ds . 

Js'(m) Jr. 

Thus, this result can also be extended to f,g £ L 2 (IR) in the sense of an L 2 (/j,)- 
limit. Hence, within the above formalism, a version of a standard Wiener process 
can be written as Wt = \l[o,t)) '/> f or t > and Wo = 0. 

To treat the occupation time of the Ornstein-Uhlenbeck process in the white 
noise framework, we need the space of Hida distributions (5)'. 

The above introduced space L 2 (/j) serves as the central space of the Gel'fand 
triple (S) C L 2 (/j) C (S)' , where (S) denotes the space of Hida test functions. 
The dual pairing of $ £ (S)' with (p £ (S) is denoted by ({ip, $)). For a detailed 
description of the construction of the Hida triple we refer to Rcf. [HKPS93 . 

Example 3.2. For a function f £ S(M), the exponential exp(i (/,•)) is an 
element of (S). 

We will characterize Hida distributions with the help of the T-transform 
and U— functionals. 

Definition 3.3 (T-transform). 

The T-transform of a Hida distribution $ £ (S)' is defined as 
T($)(/) :=«$,exp(i </,■»», 

where f £ 5(E). 

Since 1 £ (S), the expectation of a Hida distribution $ £ (S)' can be defined 

by 

E„($) :=«1,*»=T(*)(0). 

Definition 3.4 ([/-functional). 

We call F : 5(R) -> C a U -functional, if 
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1. For all f,g£ S(M), the mapping I3i4 F(g + xf) £ C is analytic and 
hence has an entire extension to C. 

2. There exist constants < K, C < oo and a continuous norm \\-\\ on S(M) 
such that 

\F(zO\<KeMc\z\ 2 U\\ 2 ) , 
for all z£<C and all £ £ S(M). 

The proof of the following equivalence theorem can be found in Rcf. [HKPS93J. 

Theorem 3.5. A mapping F : S(M) — > C is the T '-transform of a unique 
element in (S)' , if and only if F is a U -functional. 

Example 3.6. In the sense of a limit in (S)' we can define the white noise 
process as 

w(i) := (S t , oj) £ (S)' , 

where 5t denotes the Dirac delta in t > 0. This process can be considered as the 
time derivative of the Wiener process Wt(cj) = (l[ .t),w) in the sense of Hida 
distributions. 



The next result follows from Theorem 13.51 and concerns integration of a 
family of Hida distributions, see Ref. [HKPS931 IKLPSW961 IPS9T] . 

Theorem 3.7. Let (A,„4, v) be a measure space and X *— > $(A) a mapping from 
A to (S)' . We assume that the T -transform T($(A)) satisfies the following 
conditions: 

1. The mapping A i— > T($(A))(/) is measurable for all f £ S'(M). 

2. There exists a continuous norm \\-\\ on S(J&) and functions 

C £ L°° (A, v) and D £ Z/ 1 (^4, v) integrable with respect to v such that 

|T($(A))(z/)|<^(A).exp(C(A)|z| 2 ||/|| 2 ), 

for all f £ S(R), z £ C . 

Then it holds in the sense of Bochner integration in a suitable sub-Hilbert space 
of (S)' , that the integral of the familiy of Hida distributions is itself a Hida 

distribution, i.e. / $(A)c?^(A) £ (S)' and the T -transform interchanges with 
J A 

the integration 

T[ f $(A)di/(A)J = / T($(A))<fi/(A) 



Based on the above theorem, we introduce the following Hida distribution. 

Definition 3.8 (Donsker). We define Donsker's delta at x £ M corresponding 
tor] £ L 2 (R) by 

$x((v, •>) : = 7T / cx P( iX ((v, •> - x)) dX 
^7T J r, 
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in the sense of Bochner intearation \HKPS93[ \LLSW94\ [W95^ . Its T-transfo 
in f G 5(K) is given by 

T(S X ({V, '})(/) = 7= \ i exp j-r(i (V, /) - x) 2 - J (/, /) 



v/2tt (77,77) V 2(77,77) W,J/ ' 2 



Coming back to the Ornstein-Uhlenbeck process, we note that we can write 
it in the framework of white noise analysis as 

U t {u) = (U - /j) exp(-At) + n + (a exp [A(- - t)} l m (-),w) G L 2 (fi) , (3) 

for t > 0. This can be seen as follows: 

Clearly Ut is a Gaussian random variable with expectation 

E M (E/t) = / S , (R) f/ t = ^0 exp(-Ai) + /i(l - exp(Ai)) 
and covariance 

Cov(U t ,U T ) = E ll ((Ut-E M (U t ))'(U T --E ll (U T ))) 
(erexp [A(- - t)} l [0 . t) (-),o;) 

x (crcxp[A(- -r)] 1[o, t) (-),w) d/j 

2 

= ^(exp(-A-|t-r|)-exp(-A(t + r))). 

Thus, by uniqueness Ut is the Ornstein-Uhlenbeck process solving the corre- 
sponding SDE mj. 

In the special case Uq = /i = Q of the centered Ornstein-Uhlenbeck process 
= ( 7 ?t! ■)) where s %(s) = crexp [A(s — i)] lp,t)(s) € £ 2 (K), we obtain that 
the T-transform of the corresponding Donsker's delta at x G K is given by 

T(4«, t) •)))(/) = ^^oxp (-~ (f y„r )2 - \ a,/) 

v/2vr (774,774) V 2 (774,774) 2 

for / G 5(K). Using 

<r 7t ,r 7t > = . 2 jf e^-*) ds = o» 1 ~ =: * 

the expectation is readily available by 

T(^(774,-))(0) = -2=exp f-g 
V27rfc V 2fc 
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Proposition 3.9 (Expectation of occupation times). Let X t be a centered 
Ornstein-Uhlenbeck process on the time interval [0,T], where T > 0. Let 
[a, b] C M be an interval, where a = — oo and b = oo are allowed. The ex- 
pectation of the occupation time M T ^ a b ^(X t ) is given by 

E M (M TAaM {X t )) 

= I So CTf ( l-c X p(- 2 A t ) ) ~ CTf ( l-e X pT-2A t ) ) dt , (4) 

where a = \f\jo. 



Proof. The occupation time, i.e. the time the process spends in the space 
interval [a, 6] during the time T is given by 

M TAaM {X t )= [ [ 5 (X t - x) dx dt . 



Interchanging integrations due to Theorem I3.7[ we obtain the expectation of 
the occupation time of the Ornstein-Uhlenbeck process 

M M T,[a.b]( X t)) = [ [ E„(S (X t -x))dxdt 
JO J a 

T 1 f b fx 2 



exp — — J dx dt 



V^Kk J a V 2k 

If erf | —= I - erf I -= \ ,11 



, erf — h — erf — n dt, 

2 7o V y/l - exp(-2Ai) I V v/l - exp(-2At) / 



1 ' T 



with a = \f\jo~. 



4 Numerics 

4.1 Estimation of the expected occupation time by Monte- 
Carlo methods 

The expected occupation time of the Ornstein-Uhlenbeck process X t defined 
via the SDE Q can be computed using Eqn. Alternatively, one can also 
compute the occupation time using a Monte-Carlo simulation of the underlying 
process. We generate N of sample paths of the Ornstein-Uhlenbeck process 
with a fixed parameter set and compute the sample occupation time MtAo,m(') 
for each path. As in the basic idea of the Monte-Carlo simulation, the sam- 
ple average of the occupation time serves as an estimator for the expectation 
value. If large numbers of samples are considered, the estimator yields a better 
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approximation. In the sequel, we shortly outline the numerical approximation 
of a stochastic process like the Ornstein-Uhlenbeck process (HJ. 

Consider a general non-autonomous stochastic differential equation 

dX t = f(t,X t )dt + g(t,X t )dW t , X = x £R (5) 

defined in the time interval [0, T], where /, g : [0, T] x R — ► K and Wt is a 
standard Wiener process. Under mild conditions the solution of ([5]) has the 
following form 

X t = X + f f(s,X 8 )ds+ f g{s,X s )dW s , 0<t<T. (6) 
Jo Jo 

Note that the solution X t is a random variable for each t. For details on the 
existence and uniqueness of solutions to ([5]), we refer to Rcf. [Oksa07]. 

This solution can be numerically estimated by using the Euler-Maruyama 
method. We discretize the interval [0,T] using a time step At = T/L for some 
positive integer L and introduce discrete time points Tj — jAt for j — 1, 2, L. 
Let Xj denote the numerical approximation of X T . Further, we assume that the 
second integral on the right hand side of ((6]) is integrated using the Ito-version 
of stochastic integrals. Then the Euler-Maruyama method reads as 

X 3 = X j _ 1 + f( Tj ,X j _ 1 )At + gfaXj-MWr, - W Ti _ t ), ] = 1,2,..., L . (7) 

Figure Q] shows a sample path of the Ornstein-Uhlenbeck process computed us- 
ing the Euler-Maruyama method. Using the discrete version Xj of the process, 
we can easily calculate the sample occupation time My j a a(Xj). 

Remark 4.1. The accuracy of the numerical solution to the SDE can be mea- 
sured in two ways, namely strong and weak convergence. Strong convergence 
measures the accuracy on the basis of individual realizations. The weak con- 
vergence measures the accuracy of numerical methods to SDEs in case where 
the goal is to ascertain the probability distribution. For example, the Euler- 
Maruyama method has strong order of convergence 1 — \- For more details, 
see Refs. THiOTX WP~9^ . 



4.2 Direct computation of the expected occupation time 

To compute the expectation of the occupation times given by ((4]), we have to 
evaluate integrals of the type 




^g{C,\,T) (8) 
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time t 



Figure 1: A sample path of the Ornstein-Uhlenbeck process ([2]) with parameters 
A = 0.5 and a = 0.25 on the time interval [0, 25]. 



where so = C/yl - exp(-2AT). Note that s > C and s -> C for AT -> oo. 
Hence, in the limit AT — > oo, e.g. for A = 0(1) and T — »• oo, the integral gets 
singular. This together with the unbounded domain of integration poses some 
numerical difficulties, which can be overcome by splitting the integral: the part 
close to the asymptotic singularity at s = C, and intermediate part and the 
part close to infinity. We introduce s± > sq and s 2 > s± and rewrite 



g(c,x,T) 



erf (s) 



s(s 2 - C 2 ) 



ds := h 
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The second integral 1% does not pose any numerical difficulties and can easily 
be computed using Simpson's rule. However, we have to take care about the 
first and the third part. 

In l\ , we replace the error function by its quadratic Taylor polynomial 



2e" s ° 



T 2 (s) = erf(s ) H j=r- [(s - s ) - s (s - s ) 2 



at so, and get 
h 



rf(«o)C 



ds , 2e~"0 f s i (s-s )-s (s-s ) , 



1 

2C 5 



erf(so) 



2(sg + s )e- a 



1 l-C 2 /sl 
111 1-CV4 



-1 -s 2 l 

-e s o In 



(si-c)(s +C) 



s 2 -C 2 



(9) 



(10) 



Cx/i 7 x " (si+C)(s -C) 1 0F C 

The choice of si depends upon the desired accuracy of the above approximation. 



Lemma 4.2. Choosing si = s + y so(s§ — C 2 )e /or a given tolerance e > 0, 
we get an approximation error < e. 
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Proof. 

f§fi| < |erf( S )-T 2 ( S )|. r d3 

J s S{S Z - C 2 ) s <s< Sl J SQ 



s{s 2 - C 2 ) 



,d 3 erf d 3 



s eR 1 ds 3 y n 6 so{s 2 -C 2 )' 

where we introduce d = s\ — Sq and estimate the integral using the mean value 
theorem. Hence 



s ( s 2 _ (72) - 3^s (^-C 2 ) 



For the given choice of si, i.e. d = \/so(sq — C 2 ) we get 



In the third part of the integral, we replace the error-function by its limit 
erf(oo) = 1 and get 

f°° erf(s) A 1 1 g| 
L s(s 2 -C 2 ) dSRi 2C 2ln 4~C 2 ■ 

Choosing S2 > 10 yields an approximation of the error function of less than 

io- 44 . 

Figure [2] shows the expected occupation times obtained using either Eqn. Q 
or Monte-Carlo simulations. In the setting that we have shown here, both 
computational methods yield indistinguishable results. 



Remark 4.3. For large values of 2\T , the function g(C, A,T) in gets hard 
to evaluate numerically. For 2\T > 37, we obtain e~ 2XT < 10~ 16 ; the usual 
machine precision. Hence, we limit ourselves to 2XT < 37. Furthermore, in 
case of XT — ¥ oo, we obtain Sq — > C as well as Si — > C. Therefore the splitting 
of the integration introduced above does not resolve the problem with singularity 
at the bounds of the integration interval. 

4.3 Parameter estimation 

Estimating the parameters of the Ornstein-Uhlenbeck process given by ^fy based 
on Eqn. (j4]) for the occupation time is the main objective of this paper. There- 
fore, we formulate an optimization problem which we use to estimate the pa- 
rameters. 
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2 




0.2 0.4 0.6 0.8 1 

X 

Figure 2: Expected occupation times on the interval [a, b] = [—0.1,0.1], time 
horizon T = 16 and parameters a = 1 and A G (0, 1). Computations are carried 
out using either the Monte-Carlo method or Eqn. (jH). 

Let X\ t<T (t) denote the Ornstein-Uhlenbeck process with parameters A and 
a as in Eqn. Let a spatial interval [a, b) be given fixed. Then we define 
C'T,[a,b]{^, <7 ) '■= (-^r,[a.b](^A,<r(i))) as the expected occupation time of the 
Ornstein-Uhlenbeck process for time horizon T. 

Problem: Given different time horizons Ti and intervals [ctj,bj], correspond- 
ing occupation times Gi,j for i = l,...,n and j = 1, ...,m, determine the 
parameters (A, a) such that the deviation 

n m 

R(X,a) : Y,y,((h >, ':A.a;, G, f (11) 

i=i j=i 

is minimal. 

To solve this optimization problem, we apply a standard method from nu- 
merical analysis. Here, we have used the simplex search method implemented 
in Matlab as the function fminsearch, see Rcf. Mat07 . As a stopping exit 
for the optimization, we use a tolerance of 10~ 5 between successive iterations. 
To demonstrate the parameter estimation procedure, consider the following 
situation. Let T\ = 10 and Ti = 12 be the time horizons Tj for i = 1,2 
and [-0.25,0.25], [-0.5,0.5], [-0.75,0.75] [-1.0,1.0] be the intervals [a^bj] for 
j = 1,2,3,4, which use to calculate the corresponding data Gij to the above 
mentioned optimization problem. Using either the direct equation Q or the 
Monte-Carlo method, we compute occupation times Gij for the parameters 
(A, cr) = (0.15,0.90). Now we initiate the minimization procedure providing 
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the corresponding GijS for the both situations separately. The resulting esti- 
mated values are (A*,er*) = (0.150025,0.90004) in the case Gij are computed 
via direct equation fltj} and (A*,cr*) = (0.1330915,0.8764664) in the case G itj 
are computed by the Monte-Carlo methods. The following table lists more 
numerical findings that are estimated correspond to different setting. 



true parameters 


recovered from (J4| 


recovered from MC 


A 


a 


A* 


a* 


A* 


a* 


0.25 


0.75 


0.250013 


0.750024 


0.275525 


0.779474 


0.50 


0.50 


0.500012 


0.499878 


0.549382 


0.519195 


0.75 


1.25 


0.750022 


1.249987 


0.720177 


1.225472 


1.00 


2.00 


0.999917 


2.000018 


1.026446 


2.026084 


1.25 


2.50 


1.250011 


2.500011 


1.309042 


2.554252 



The parameters recoverd from occupation times generated using £1|) (columns 3 
and 4) agree better than those recoverd from occupations times generated with 
the help of Monte-Carlo simulations (last two columns). This is not surpris- 
ing, since we used the same underlying equation to generate and to recover the 
parameters. However, also for the parameters covered from the Monte-Carlo 
simulations we have a difference of about 10% between the true and the re- 
covered data. This is for most applications a sufficient accuracy. Nevertheless, 
increasing the number of samples in Monte-Carlo simulation we can improve 
the accuracy of the recomputed parameters. 

5 Conclusion 

We derived an equation to compute directly the expected occupation time of the 
centered Ornstein-Uhlenbeck process. This allows to identify the parameters of 
the Ornstein-Uhlenbeck process for available occupation times via a standard 
least squares minimization. To test our method, we generated occupation times 
and recovered the parameters with the above mentioned procedure. Within the 
range of our numerical experiments, we found very good agreement. This gives 
hope to be able to estimate parameters in industrial fleece production processes 
from measurable quantities like the mass per area. However, to get closer to the 
industrial applications we have to extend our method to the 2d-case and more 
involved processes than the standard Ornstein-Uhlenbeck process. 
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